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1 Introduction 

Reliable estimation of the local spatial extent of 
neural activity is a key to the quantitative analysis of 
MEG sources across subjects and conditions. In 
association with an understanding of the temporal 
dynamics among multiple areas, this would 
represent a major advance in electrophysiological 
source imaging. 

Parametric current dipole approaches to MEG (and 
EEG) source localization can rapidly generate a 
physical model of neural current generators using a 
limited number of parameters. However, 
physiological interpretation of these models is often 
difficult, especially in terms of the spatial extent of 
the true cortical activity. In new approaches using 
multipolar source models [3, 5], similar problems 
remain in the analysis of the higher-order source 
moments as parameters of cortical extent. Image- 
based approaches to the inverse problem provide a 
direct estimate of cortical current generators, but 
computationally expensive nonlinear methods are 
required to produce focal sources [1,4]. Recent 
efforts describe how a cortical patch can be grown 
until a best fit to the data is reached in the least- 
squares sense [6], but computational considerations 
necessitate that the growth be seeded in predefined 
regions of interest. 

In a previous study [2], a source obtained using a 
parametric model was remapped onto the cortex by 
growing a patch of cortical dipoles in the vicinity of 
the parametric source until the forward MEG or 
EEG fields of the parametric and cortical sources 
matched. The source models were dipoles and first- 
order multipoles. 

We propose to combine the parametric and imaging 
methods for MEG source characterization to take 
advantage of (i) the parsimonious and 
computationally efficient nature of parametric 
source localization methods and (ii) the anatomical 
and physiological consistency of imaging techniques 
that use relevant a priori information. By 
performing the cortical remapping/imaging step by 
matching the multipole expansions of the original 


parametric source and the equivalent cortical patch, 
rather than their forward fields, we achieve 
significant reductions in computational complexity. 

2 Methods 

2.1 Parametric source imaging 

Though the equivalent cortical patch may be seeded 
with the solutions from any parametric localization 
method, we routinely use Recursively Applied and 
Projected (RAP)-MUSIC which does not require 
that the number of sources be known a priori, as is 
the case for least squares fitting [7], The version 
available in Brainstorm (http://neuroimage.usc.edu) 
estimates parameters for sources up to quadrupolar 
complexity. 

2.2 Seeding the cortical growth 

The first step consists of finding one or several best 
seeding locations from which to start growing a 
patch. We select the initial cortical locations for 
region growing using areas in which the local 
current density is larger than a predefined threshold 
in a Minimum-Norm (MN) image of the cortical 
currents. 

Let G 0 be the gain vector characterizing the 
estimated source J o at location r o for which an 
equivalent cortical patch is to be found. The cortical 
surface from a subject MRI is tesselated and an 
elemental current dipole allowed at each vertex on 
this surface. A gain vector G is assigned to the 
dipole at location r ; oriented normally to the 
cortical surface. Finding a cortical equivalent to J 
consists of finding the set of elemental cortical 
dipoles £ and their amplitudes y = [y,, i e £] so 
that G 0 =G ? y + b = ^G,.y l . +b, where b is the 

error due to the discrepancy between the parametric 
source and its cortical counterpart. £ is intialized by 
thresholding a regularized MN estimate of y, 



y = argmin |G 0 - G ? y ||“ + A||y|| 2 : 

£ = {ie C, |j>,-| >T }. (1) 

where C denotes the entire set of N cortical vertices. 

2.3 Iterative estimation of the source 
extension 

The patch growth consists of a recursive estimation 
of the local current density on a growing number of 
source candidates in the vicinity of every seed in £ 
until a functional U[ is minimized (for seed i at 
iteration k). The functional U\ is defined in the 
following section. 

A. For every seed z'e £, 

1. Initialization: set k= 1 and the patch around 
source i to ■&[ = {/} and U' 0 = 0; 

2. Estimate y,., then compute U ‘ k ; 

3.If | U[ !>£[/'_, , 

(i) Grow the patch by including the vertices 
connected to the source(s) in , 

(ii) Set k=k +1 and go back to A.2; 
else: 

(i) Define U‘ = U [_ x ; 

(ii) Define the best patch obtained from seed i, 

rT=#'_i ; 

(iii) Proceed to the next seed. 

B. Patch Merging: 

Select a subset of patches {lT,ze /} with 
l = {i,U‘ <{(u J )-3o(u '))}, where (u 1 ) (res. 
a ([/ 1 )) is the mean (res. the standard deviation) of 
the values of the functionals in e £ }. The 

final patch II is then defined as: n = |J IT • 

For every seed i and at every iteration step k, the 
source amplitudes - which may be non-null only 
within ■&[ - are estimated using a Bayesian estimator 
[1][4] as described in the following. 

2.4 Bayesian imaging of the local current 
density 

The local current density is modeled as a random 
field using extensions of models described in [4], 
We characterized the current density at every vertex 
with a continuous random variable of dipole 
amplitude z, and a binary indicator process x,- to 
model whether source i is on or off. The local 
equivalent source density can thus be written as 

y. = x. * Z.. 


A MAP Bayesian estimator of y is given by: 

- - „ p{b\x,z)p(x)p(z) 

y = Xz = argmax- 1 -^-, (2) 

Mb) 

where X = diag(x), and x and z are assumed to be 
independent. We assume that the source amplitudes 
z follow a zero-mean Gaussian process and that the 
indicator process x follows the assumption that 
cortical activity is focal and sparse. We hence define 
p{\) to be a Gibbs distribution with the associated 
energy v(\) function given by: 


K (x)=IJ a j x j + PjY,~ 


Yj U 


( 3 ) 


where the parameters a >0 and /? > 0 determine 
the relative weights of the sparseness and clustering 
terms; v j is the set of nearest neighbors of vertex j; 
and Yj„ > 0 is a coefficient translating the anatomic- 


functional priors between locations j and u. y ju is 
proportional to the distance between the 2 vertices 
and to the discrepancy between their orientations 
(see [8] for details). 

The global functional to minimize for seed i at 
iteration k is thus: 

D'(x,z) = ^-[G 0 -G ? Xz] t C“‘[G 0 -G f Xz] 

+ ~z T C) 1 z + F(x) 

( 4 ) 

where C n (res. C z ) is the covariance matrix for the 
perturbation process b (res. the amplitude process z). 
U‘ k (x,z) is minimized using an approach based on 
the technique described in [4] with an ICM (iterated 
conditional mode) optimization of the binary 
indicator process. 


2.5 Equivalent cortical patch estimation 
using scalar multipole expansions 

A direct approach to remapping inspired by inverse 
problem techniques would be to define the gain 
vector G as the forward field of source i sampled at 
the sensor array. But there is no fundamental 
requirement that we directly match the sampled 
magnetic field B of the equivalent parametric and 
cortical sources. Indeed, as the current distribution is 
zero outside the head, B is the negative gradient of a 
scalar magnetic potential V m obeying Laplace’s 
equation. V m can hence be expressed in terms of a 
spherical harmonic multipole expansion, the 
moments of which are related to the current source 
density. Hence matching the magnetic fields 



everywhere outside the head reduces to matching the 
moments of the original source density. If the 
original source consists of dipoles for instance, it is 
straightforward to relate it to a multipole series 
about the origin [9], Shift equations also exist to 
relate a multipolar expansion about a given point to 
one expanded also about the origin [10]. As 
potentials related to moments of order n fall off as 
r“ ( " +1) , we may assume that most of brain sources 
may be reasonably described with moments of order 
up to 4 (i.e. octapoles). This means that the cortical 
remapping can be done with gain vectors describing 
only the first 15 moments of the scalar multipole 
expansion of the source about the origin (dipole: 3; 
quadrupole:5; octapole:7). Hence for every original 
and cortical source, we compute the following gain 
vector containing the 15 first moments of the 
multipole expansion of the source about the origin: 



where m = ^ r x J (r) is the equivalent magnetic 
dipole of source j(r) located at r; Q'are the 
quadrupole moments of the source, and O' its 
octapolar moments, both expanded about r; A ; arc 
matrix kernels derived for the shift equations (see 
[10] for details); and I, stands for the identity matrix 
of size j. The matrix on the left in (5) - subsequently 

denoted by -^1. - allows a scaling of all the 

multipole moments to the same unit and facilitates 
better conditioning of G ? . For dipoles, the 
multipole expansion is linear in terms of the 
equivalent magnetic moment and simplifies to: 

G=^I. A, [rxo,], (6) 

r A 2,l 

o, being the orientation of source i. 

The advantages of this methodology include: (i) 
original and cortical sources are explicitly matched 
without reference to the source fields; (ii) there is no 
need to compute the forward fields on the entire set 
of several thousand cortical vertices, thus saving 
computation time, memory and storage; (iii) the 
remapping and imaging approaches also benefit 
from faster computational speed as the gain matrices 
to manipulate are 15 rows tall rather than having 


dimension equal to the number of sensors, which 
could be several hundred for recent MEG systems 
(iv) patches need not have uniform activation fields. 
Our Matlab code for remapping on a 5,000-node 
cortex takes about 5secs on a PHI 650MHz, 128MB. 

3 Results 

3.1 Simulations 

A 5,000-node white/gray matter interface was 
extracted using the BrainSuite software (USC, 
http://neuroimage.usc.edu). Monte-Carlo (MC) 
simulations were performed by growing about 2300 
cortical patches at randomly selected locations on 
the cortical surface with areas ranging from 0.3cm 2 
to 15.3cm 2 (median 4.1cm 2 ). Uniform illumination 
was assigned to the cortical dipoles within a patch 
using a 100-time-sample waveform for active 
dipoles. MEG signals were simulated on 138 axial 
gradiometers (5cm baseline) uniformly distributed 
about the upper hemisphere of a spherical head. 
Gaussian white noise was added to the signals with a 
uniform level across all the channels of 10% of the 
peak of maximum amplitude. 

At every MC trial, the parametric estimation with 
RAP-MUSIC consisted of finding the best dipole 
source that could account for the signal generated by 
the patch. Cortical remapping was then applied to 
the source found. C n (res. C z ) is chosen as <r 2 I (res. 
ol \) with a 2 = 10' 2 (res. a 2 = 10[nA.m] 2 ). a. and 
were set to 10' 5 for every source, and no priors 
besides connectivity were taken into account when 
generating the original patches, hence y ] U =1 for all 
pairs of neighbors. £ was set to 10“ 6 . 

3.2 Results 

The patches generated in the MC simulations were 
gathered in 5 classes according to their area. Each 
class was labeled by the average value of the patch 
areas within that class. Following source remapping, 
the estimated area and the subspace correlation with 
the original cortical patch were computed. We 
designed an area-matching criterion from the 
original and the remapped patches: it consisted of a 
uniformly weighted sum of (i) the relative 
localization error of the patch centroid; (ii) the 
discrepancy with the original patch area; (iii) the 
overlap with the original patch and (iv) an affine 
transform of the resulting subspace correlation sc: 
(10 sc — 9)e [0,1]. This area-matching criterion takes 
its values between 0 (no match) and 100% (perfect 
match). Fig.l shows there is no significant 



degradation of the remapping with increasing area of 
the original cortical patch (average matching 
criterion: 87%). 



Figure 1: Area-matching criterion for each of the 5 
patch classes. 

More specifically, the subspace correlation with the 
original patch was excellent (average 0.993) even 
though the original parametric source tended not to 
fit the original forward field as well when the active 
cortical area increases (average 0.985) (see Fig.2). 



Figure 2: Comparison of subspace correlation 
between the original patch and the parametric 
source (squares) or the remapped source (circles). 

A nice feature of the remapping method is that the 
original area is restored with surprising accuracy as 
shown Fig. 3. 



Figure 3: Estimated (circles) vs. true (squares) value 
of the surface of the active cortical surface. 

The estimator succeeds in recovering quantitatively 
the area of the original surface with a correlation of 
0.99. 

4. Conclusion 

We have presented a fast and promising method to 
localize and estimate the spatial extent of activated 
cortical regions in MEG. Further developments will 
consist of evaluating and correlating these new 
quantitative indices under various protocols and 
conditions. 
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